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Executive  Summary 


A  Presidential  Mandate  in  2013  committed  the  Department  of  Defense  (DoD)  to 
converting  at  least  20%  of  its  energy  demand  over  to  renewable  energy  sources  by 
2020,  to  include  solar,  wind,  water,  and  geothermal.  The  Army  also  declared  its 
intent  on  becoming  a  net-zero  energy  consumer  by  2030.  In  order  to  accomplish 
these  goals,  renewable  energy  sources  need  to  be  integrated  into  existing  electrical 
power  control  architecture.  One  method  is  to  integrate  renewable  energy  resources 
into  the  traditional  tactical  microgrids,  creating  hybrid  microgrids.  The  challenge 
with  such  a  task  is  the  inconsistency  in  the  power  generated  by  renewable  resources. 

In  order  to  create  an  effective  and  functioning  hybrid  microgrid,  control  algorithm 
is  needed  to  regulate  the  amount  of  power  coming  into  and  going  out  of  power  grid 
components.  For  this  high-level  controller  to  work  properly,  it  needs  to  quantify  the 
amount  of  power  produced  by  the  renewable  resource  that  will  be  available  at  the 
site,  before  the  power  is  drawn  off  to  its  designated  load.  This  problem  can  be 
abated  with  the  implementation  of  an  atmospheric  model,  such  as  a  Solar  Radiation 
Flux  (SRF)  Model.  This  report  documents  a  program  coded  for  renewable  energy 
applications  that  was  based  on  Ralph  Shapiro’s  SRF  Model,  which  was  selected 
based  on  its  use  of  insitu-only  measurement  input. 

There  are  numerous  factors  that  limit  the  total  amount  of  solar  radiation  flux 
received  at  the  ground  level  after  traversing  the  atmosphere.  The  model 
documented  in  this  report  simplifies  these  variables  and  implements  them  in  an 
algorithm  that  predicts  the  net  solar  radiation  flux  received  at  ground  level 
throughout  the  entire  solar  day.  The  model  represents  the  atmosphere  as  3 
homogeneous  layers,  and  quantifies  the  effects  on  the  vertical  radiation  in  terms  of 
reflectance,  transmittance,  and  absorption.  These  3  terms  are  a  function  of  cloud 
type,  thickness,  and  current  atmospheric  conditions. 

The  model’s  purpose,  necessity,  and  development  are  detailed  in  Chapter  1. 
Chapter  2  describes  the  model  equations  and  data  tables,  which  are  summarized  in 
a  model  algorithm  flowchart.  A  preliminary  comparison  between  the  model  output 
and  measured  data  sampled  by  Atmospheric  Renewable  Energy  Field  Study  #2 
(ARE2)  pyranometers  is  provided  in  Chapter  3.  Chapter  4  summarizes  model 
results  and  provides  recommendations  for  future  tactical  energy  unit  independence 
applications. 

One  recommendation  for  future  model  development  is  to  create  a  stepwise  function 
to  facilitate  the  insertion  of  insitu  meteorological  observations  throughout  a  single 
day.  These  perpetual  and  timely  solar  irradiance  updates  could  yield  a  more 
accurate,  dynamic  model.  With  such  a  model,  the  isolated  power  grid  controllers 
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would  be  well  poised  for  exploiting  and  optimizing  the  atmospheric  contributions 
to  a  hybrid  power  distribution. 
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1.  Introduction 


In  2013,  a  Presidential  Mandate  committed  the  DoD  to  convert  at  least  20%  of  its 
energy  demand  over  to  renewable  energy  sources  by  2020,  to  include:  solar,  wind, 
water,  and  geothermal.  Additionally,  the  Army  is  intent  on  becoming  a  net-zero 
energy  consumer  by  2030.  In  order  to  accomplish  this,  the  Army  needs  to  find  a 
way  to  implement  renewable  energy  sources  into  existing  power  control 
architecture.  Solar  generated  power  is  one  renewable  energy  resource  that  could  be 
used  to  accomplish  the  Army’s  goals. 

Integrating  solar  power  into  traditional  isolated  microgrids’  is  challenged  by  the 
variability  of  solar  radiation  fueling  the  power.  While  designing  hardware 
alternatives  into  the  grid  configuration  can  counter  this  issue,  for  a  relatively 
constant  electrical  power  source  one  can  also  develop  a  model  that  will 
anticipate/predict  the  total  amount  of  surface  solar  radiation  flux.  This  solar 
radiation  flux  can  then  be  combined  with  the  total  surface  area  of  the  solar  energy 
sources,  namely  photovoltaic  (PV)  panels,  to  roughly  determine  the  energy 
producing  potential  of  an  installation’s  solar  array.  The  implicit  advantage  of  using 
a  predictive  model  is  that  it  would  also  allow  for  an  estimate  of  total  energy  that 
could  be  collected  and  stored  over  the  entire  day  in  order  to  more  effectively  ration 
the  fossil  fuel  sources  that  the  installation  has  available.  The  “ideal”  model  would 
be  location-specific.  Thus,  it  would  be  able  to  support  utility/installation,  microgrid, 
and  tactical  scale  power  resources. 

1.1  The  Challenge  and  Solution 

The  challenge  with  implementing  renewable  energy  sources  into  traditional 
microgrids,  making  them  “hybrid  microgrids,”  is  the  inconsistency  in  the  amount 
of  power  generated  by  the  renewable  resource.  In  order  to  create  an  effective  and 
functioning  hybrid  microgrid,  there  needs  to  exist  some  form  of  control  algorithm 
that  regulates  the  amount  of  power  coming  into  and  going  out  of  the  power  grid 
components.  For  this  high-level  controller  to  work  properly,  it  needs  to  quantify  the 
amount  of  power  produced  by  the  renewable  resource,  which  for  this  research  is 
solar,  that  will  be  available  at  the  site  before  the  power  is  drawn  off  to  its  designated 
load.  This  problem  can  be  abated  with  the  implementation  of  a  model,  such  as  a 
Solar  Radiation  Flux  (SRF)  Model.  The  SRF  Model  has  the  potential  of  providing 
atmospheric  input  that  grid  operators  can  use  to  accurately  predict  the  net  solar 
radiation  flux  received  throughout  the  solar  day  after  taking  routine  surface 

Microgrid  is  defined  as  “multiple  power  resources  assembled  as  a  single  system  (generator, 
storage,  distribution  and  load),  with  the  ability  to  run  independently  as  an  “island”  and/or,  as  an 
integrated  part  of  a  larger  grid.”(Vaucher  et  al  2016) 
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meteorological  observations,  location  information,  and  time  infonnation  into 
account.  The  model  chosen  to  answer  the  challenge  was  primarily  based  on  work 
published  by  Ralph  Shapiro  (1982).  The  Shapiro  SRF  Model  was  selected  due  to 
its  ability  to  use  insitu-only  measurements  input.  This  report  documents  the  results 
of  translating  Shapiro’s  written  concepts  into  a  computer  code,  along  with  several 
improvements  to  strengthen  and  retailor  the  code,  in  preparation  for  its  renewable 
energy  applications. 

1.2  The  Solar  Radiation  Flux  Model 

When  one  considers  solar  power  as  a  relatively  constant  source,  as  the  Sun  will  rise 
and  set  every  day  in  most  locations  around  the  world,  the  concern  then  becomes 
quantifying  the  various  effects  that  filter  solar  radiation  flux  before  it  reaches  the 
ground.  The  model  encoded  represents  these  factors  as  the  relative  reflectance, 
transmittance,  and  absorption  of  atmospheric  layers  that  will  act  on  the  solar 
radiation  as  it  traverses  strata.  These  tenns  are  a  function  of  cloud  type,  size,  and 
density. 

To  create  a  practical  computational  environment  in  which  enough  data  can  be 
collected  on  each  of  the  model’s  key  tenns,  some  aspects  of  the  atmosphere  and 
clouds  must  be  standardized.  Theoretically,  the  atmosphere  can  be  divided  into  any 
number  of  descriptive  layers  with  slightly  different  qualities  in  which  different 
cloud  types  would  reside.  For  the  sake  of  simplicity,  this  model  organizes  the 
atmosphere  into  3  homogeneous  layers  to  mirror  the  general  vertical  atmospheric 
structure  used  by  trained  Sky  Observers  or  Aerographers  (Orvis  et  al.  1984).  With 
the  wide  variation  of  cloud  types  and  cloud  microphysics,  clouds  are  classified  into 
3  groups  to  populate  the  SRF  Model:  top-layer  clouds  (thick  or  thin  cirrus, 
cinostratus,  cirrocumulus),  mid-layer  clouds  (altostratus  and  altocumulus),  and 
low-layer  clouds  (cumulus,  cumulonimbus,  and  stratocumulus  or  stratus). 

This  infonnation,  when  combined  with  the  relative  position  of  the  Earth  and  Sun 
detennined  from  site  location  and  date,  can  be  used  to  find  the  net  solar  radiation 
flux  received  at  ground  level.  When  propagated  over  an  entire  solar  day,  the 
inadiance  takes  on  a  Gaussian  curve  distribution  that  starts  at  the  local  sunrise  and 
ends  at  the  local  sunset.  The  relatively  stable  results  of  the  model  appeared  to  match 
pyranometer  solar  radiation  data  collected  under  clear  skies  during  the  Atmospheric 
Renewable  Energy  Field  Study  #2  (ARE)  at  White  Sands  Missile  Range,  NM. 
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2.  Procedure 


In  order  to  develop  a  programmable  model,  it  was  necessary  to  organize  the  model 
into  components  that  could  be  coded  into  MATLAB:  inputs,  constants,  functions, 
and  outputs,  as  reflected  in  Fig.  1 .  After  creating  the  algorithm  flowchart  based  on 
Shapiro  (1982),  it  became  clear  that  the  final  output,  solar  radiation  flux  (X3),  was 
a  product  of  2  separable  elements — the  extraterrestrial  radiation  (Xo)  and  the 
coefficient  of  transmissivity  due  to  atmospheric  conditions  (X3).  Extraterrestrial 
radiation  is  represented  in  Fig.  1  by  the  elements  below  the  red  line,  and  the 
coefficient  of  transmissivity  is  represented  by  the  elements  above  the  red  line. 
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longitude 


Fig.  1  Algorithm  flowchart 


2.1  Extraterrestrial  Radiation 


Extraterrestrial  radiation,  measured  in  W/m2,  is  a  measure  of  the  solar  radiation 
received  at  the  top  of  the  atmosphere  represented  by  Eq.  1.  (Shapiro  1982) 

X0  =  S  cos(z)  (1) 

fd\2 

Where  S  is  the  quiet  Sun  solar  constant,  (  -J  is  a  function  of  the  ellipticity  of  the 

Earth’s  orbit  and  the  position  of  the  Earth  in  its  orbit  around  the  Sun  represented  by 
Eq.  2,  and  z  is  the  zenith  angle  represented  by  Eq.  3.  (Shapiro  1982)  During 
minimum  solar  activity,  S  is  roughly  1369.2  W/m2. 

(f)2  =  [1.000140  +  0.016726 cos  (2) 

JD  is  the  local  Julian  Date,  the  number  of  days  since  January  1  to  the  current  date, 
and  day  is  the  number  of  days  in  a  year,  which  changes  from  365  to  366  on  leap 
years. 

z  =  acos(sin(9)sin(8 )  +  cos(9)cos(8)cos(h))  (3) 

0  is  the  current  latitude,  6  is  the  declination  angle  described  in  Section  2.1.1,  and  h 
is  the  hour  angle  described  in  Section  2.1.2. 


2.1.1  Declination  Angle 


Declination  angle  is  the  seasonal  tilt  of  the  Earth  on  its  axis.  This  tilt  affects  the 
total  amount  of  solar  radiation  received  by  certain  locations.  The  process  to 
calculate  the  declination  angle  is  outlined  in  Eqs.4-13.  (Meeus  1998) 


JC 


(JD  -2451545) 
36525 


GML 


280.46646  +/C+(36000.76983  +  0.0003032/C) 
360 


) 


GMA  =  357.52911  +  /C(35999.05029  -  0.0001537  *  JC) 


EEO  =  0.016708634  -/C(0.000042037  +  0.0000001267/C) 


(4) 

(5) 

(6) 
(7) 


EOC  =  sin(GMA) (1.914602  -/C(0.004817  +  0.000014/C  + 
sin(2GMA)(0.019993  -  0.000101/C)  +  0.000289sin(3CM71)))  (8) 

TL  =  GML  +  EOC  (9) 


AppL  =  TL  —  0.00569  -  0.00478sin(125.04  -  1934.136/C)  (10) 
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MObl  =  23  + 


(ID 


26+  21.448- 


7C(46.815+yC(0.00059-0.001813yC))N 


60 


OblCorr  =  MObl  +  0.00256cos(125.04  -  1934.136 JC)  (12) 

8  —  asin(sin(OblCorr)sin(AppL ))  (13) 

The  absolute  Julian  Day  (JD )  is  a  count  of  the  number  of  days  since  noon  on 
January  1,  1900  on  the  Gregorian  calendar.  Absolute  Julian  Day  is  converted  into 
Julian  Century  (JC).  The  Geometric  Mean  Longitude  of  the  Sun  (GML)  and 
Geometric  Mean  Anomaly  (GMA)  are  used  to  calculate  the  position  of  a  body  in 
an  elliptical  orbit  as  a  function  of  JC.  Eccentric  Earth  Orbit  (EEO)  is  a  measure  of 
the  amount  that  Earth  deviates  from  a  perfect  circle  in  its  orbit  about  the  Sun. 
Equation  of  Center  of  the  Sun  (EOC)  is  the  angular  difference  between  the  position 
of  a  body  in  its  elliptical  orbit  and  the  position  it  would  occupy  if  motion  were 
uniform.  Sun's  True  Longitude  (TL)  is  the  ecliptic  longitude  at  which  an  orbiting 
body  could  actually  be  found  if  its  inclination  were  zero.  Sun's  Apparent  Longitude 
(AppL)  is  corrected  for  apparent  displacement  of  a  celestial  body  based  on  the 
Earth’s  moving  frame  and  the  rocking  motion  that  takes  place  about  Earth’s  axis  of 
rotation,  as  opposed  to  mean  longitude.  MObl  and  oblique  correction  (OblCorr) 
account  for  the  angle  at  which  the  Earth  is  tilted,  toward  or  away  from  the  Sun 
depending  on  date,  about  its  poles.  (Meeus  1998) 


2.1.2  Hour  Angle 

Hour  angle  (h)  is  a  measure  of  the  Sun’s  angle  to  the  position  of  the  observer 
starting  at  solar  noon,  0°,  and  progressing  in  both  directions  by  15°  per  hour  for  the 
entire  solar  day.  The  angles  are  negative  to  the  East,  toward  sunrise,  and  are  positive 
to  the  West,  toward  sunset. 


Y  = 


(  OblCorr \ 

tan  - 

V  2  ) 


(14) 


EOT  =  4(Ysin{2GML )  -  2EEOsin(GMA )  + 
4EEO(Y)sin(GMA)cos(2GML )  -  O.SY2sin(4GML)  - 
1.25EE02  sin(2GMA)) 

(15) 


SRAng  =  acos 


(  cos( 90.833) 
\  cos(9 ) cos(S ) 


—  tan(9)tan(8)  J 


(16) 


SN  = 


(720— 4long— EOT +60tz) 
1440 


(17) 


SR  =  SN 


(18) 
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1440 


(19) 


SS  =  SN  + 


4 SRang 
1440 


SL  =  8  SR 


TST  =  mod 


(  0.l(l440) 
24 


+EOT+4long-60tz-60 


1440 


h 


TST 

=  —  +180 
4  — 


(20) 

(21) 

(22) 


The  term  Y  is  a  placeholder  for  the  modification  of  the  oblique  correction 
(OblCorr),  which  puts  it  in  a  form  usable  by  equations  going  forward.  The  Equation 
of  Time  (EOT)  is  an  adjustment  to  the  clock  based  on  the  orbit  of  the  Earth  about 
the  Sun  during  a  single  year.  This  term  modifies  the  time  of  Solar  Noon  based  on 
the  JD.  Sunrise  Angle  (SRAng)  is  the  angle  at  which  the  Sun  rises  above  the 
horizon.  Solar  Noon  (SN)  is  when  the  Sun  crosses  the  local  standard  meridian  and 
is  at  the  highest  point  in  the  sky.  Sunrise  (SR)  and  Sunset  (SS)  are  the  adjust  times 
at  which  the  Sun  crosses  the  Eastern  and  Western  horizons,  respectively.  Sunlight 
Length  (SL)  is  the  total  duration  or  time  the  Sun  is  above  the  horizon.  True  Solar 
time  (TST)  is  the  time  based  on  the  apparent  solar  day.  In  the  calculation  of  h,  when 
TST  is  less  than  zero,  then  180  is  added. 


2.2  Reflectivity,  Transmissivity,  and  Absorption 

Reflectivity,  transmissivity,  and  absorption  are  the  tenns  used  to  modify  the 
extraterrestrial  radiation  into  the  amount  actually  received  at  ground  level.  The 
model  code  aligns  with  Shapiro’s  approach,  by  presuming  that 

Ra  +  Ta  +  Aa  =  1 ,  (23) 

Where  R  is  reflectivity,  T  is  transmissivity,  A  is  absorption,  and  k  is  representative 
of  the  layer  that  is  being  calculated.  A  simple  2-stream  approximation  was  then 
used,  assuming  that  the  radiative  transport  of  solar  flux  only  travels  toward  and 
away  from  the  surface  without  a  major  scattering  effect.  Shapiro  calculated  mean 
climatological  reflection  and  transmission  coefficients  using  the  National  Climatic 
Data  Center’s  hourly  solar  radiation  and  meteorological  observations  SOLMET 
data  resource.  His  results  describe  the  overall  transmissivity  and  reflectivity  of  each 
atmospheric  layer.  The  processed  solar  radiation  assumptions  are  described  using 
Eqs.  24  and  25,  respectively  (Shapiro  1982;  SOLMET  1978).  The  equations  serve 
to  create  a  decimal  value  that  represents  the  amount  of  solar  radiation  reflected  or 
transmitted  based  on  current  atmospheric  conditions. 


Rk  =  <PkPk  +  (1  -  (pk)rk 

(24) 

Tk  =  (Pk^k  +  (1  -  (Pk)tk 

(25) 
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Where  (p  is  the  a  weighting  coefficient  expanded  upon  in  Eq.  26  and  27,  p  is  the 
overcast  layer  reflectivity  coefficient  calculated  from  Eq.  27,  r  is  the  clear  layer 
reflectivity  coefficient  calculated  by  Eq.  27,  x  is  the  overcast  layer  transmissivity 
coefficient  calculated  from  Eq.  28,  and  t  is  the  clear  sky  transmissivity  coefficient 
calculated  from  Eq.  28. 


<Pk  =  Wfk  (26) 

W  —  C0  +  C±cosz  +  C2fk  +  C3fkcosz  +  C4  cos2  z  +  C5f2K  (27) 

The  weighting  coefficient  (cpk)  is  a  function  of  the  fractional  cloud  coverage  (fk)  per 
layer,  a  biquadratic  polynomial  function  of  weight  (W)  for  each  cloud  type  and 
fractional  cloud  coverage,  and  the  solar  zenith  angle  (z).  The  C  series  coefficients, 
shown  in  Table  1,  are  a  result  of  Shapiro’s  observations  and  provide  an  adjusted 
measurement  for  each  cloud  type.  In  this  way,  both  the  overcast  and  clear 
conditions  for  each  layer  are  taken  into  consideration  in  Eqs.  23  and  24  to  determine 
the  overall  reflectance  and  transmission. 


Table  1  C-series  coefficients  (Shapiro  1982) 


Co 

Ci 

C2 

Cs 

C4 

Cs 

Thin  Ci/Cs 

0.675 

-3.432 

1.929 

.842 

2.693 

-1.354 

Thick  Ci/Cs 

1.552 

-1.957 

-1.762 

2.067 

.448 

.932 

As/Ac 

1.429 

-1.207 

-2.008 

.853 

.324 

1.582 

Sc/St 

.858 

-1.075 

-.536 

.750 

.322 

.501 

Cu/Cb 

2.165 

-1.277 

-3.785 

2.089 

-.387 

2.342 

Overcast  and  clear  layer  reflectivities  are  calculated  with  the  same  equation,  Eq. 
28,  using  “a”  series  coefficients  (see  Table  2)  that  are  the  result  of  Shapiro’s 
observations.  Here  the  F/K  condition  represents  when  there  is  fog  or  smoke  present. 

rk  or  pk  =  a0  +  aqcos^  +  a2  cos2  £  +  a3  cos 3  £  (28) 


Table  2  “a”  series  coefficients  (Shapiro  1982) 


Rk 

ao 

ai 

a2 

a3 

n 

.12395 

-.34765 

.39478 

-.14627 

1'2 

.15325 

-.39620 

.42095 

-.14200 

T3 

.15946 

-.42185 

.48800 

-.18493 

r3  (F/K) 

.27436 

-.43132 

.26920 

-.00447 

pi  (Thin) 

.25674 

-.18077 

-.21961 

.25272 

pi  (Thick) 

.60540 

-.55142 

-.23389 

.43648 

p2  (As/Ac) 

.66152 

-.14863 

-.08193 

.13442 

p3  (Sc/ St) 

.67072 

-.13805 

-.10895 

.09460 

p3  (Cu/Cb) 

.71214 

-.15033 

.00696 

.03904 
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Overcast  and  clear  layer  transmissivities  are  calculated  with,  Eq.  29,  using  different 
“b”  series  coefficients  that  are  the  result  of  Shapiro’s  observations  (see  Table  3). 

tk  or  rk  —  b0  +  fqcos^  +  b2  cos2  £  +  b3  cos3  £  .  (29) 

Table  3  “b”  series  coefficients  (Shapiro  1982) 


Tk 

bo 

bi 

b2 

b3 

ti 

.76977 

.49407 

-.44647 

.11558 

t2 

.69318 

.68227 

-.64289 

.17910 

t3 

.68679 

.71012 

-.71463 

.22339 

t3  (F/K) 

.55336 

.61511 

-.29816 

-.06663 

ii  (Thin) 

.63547 

.35229 

.08709 

-.22902 

ii  (Thick) 

.26498 

.66829 

.24228 

-.49357 

12  (As/Ac) 

.19085 

.32817 

-.08613 

-.08197 

13  (Sc/St) 

.17960 

.34855 

-.14041 

.00952 

13  (Cu/Cb) 

.13610 

.29964 

-.14875 

.01962 

There  is  one  additional  case  in  which  the  values  of  r,  p,  t,  and  x  are  taken  from  a 
table,  and  not  calculated  by  their  respective  equations.  This  special  case  occurs 
when  the  layer(s)  above  the  layer  being  calculated  are  considered  overcast,  fk-i  > 
9/10ths.  In  this  instance,  the  current  layer  is  considered  to  be  diffused,  meaning  that 
all  radiative  fluxes  are  no  longer  considered  direct.  As  such,  the  relative  reflectivity 
and  transmissivity  are  taken  from  Table  4  in  order  to  coincide  with  the  difference 
in  the  amount  of  radiation  allowed  to  enter  each  layer.  The  values  in  Table  4  are  a 
result  of  Shapiro’s  observations  on  the  basis  of  increased  path  length  as  compared 
to  a  unit  path  length  of  the  cosine  of  the  zenith  angle. 

Table  4  Diffused  coefficients  (Shapiro  1982) 


Clear  layer  (ru,  tk) 

Overcast  layer  (pk,Tk) 

Ri 

r2 

.040 

.560 

r3 

.045 

.520  (Cu/Cb) 

r3 

.116  (F/K) 

.609  (Sc/St) 

Ti 

t2 

.905 

.361 

t3 

.900 

.400  (Cu/Cb) 

t3 

.788  (F/K) 

.311  (Sc/St) 

Finally,  the  coefficient  of  transmissivity  due  to  atmospheric  conditions  is  calculated 
by  Eq.  30,  and  multiplied  with  extraterrestrial  radiation  (Xo)  to  yield  solar  radiation 
flux  at  the  ground  level  (X3). 
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*3  =  T1T2T3/D2  (30) 

D2  =  (1  -  R2Rg)[a  -  R1R2X 1  -  «2«3)  -  R1R3T2)]  -  RgTitt  1  -  + 

R\T2]  (31) 

In  Eq.  3 1 ,  Rg  represents  the  relative  albedo  or  reflectivity  of  the  ground  and  foliage 
layer  represented  by  a  decimal. 

A  similar  SRF  Model  developed  by  David  Tofsted  calculated  the  surface  solar 
radiation  flux  for  individual  times  throughout  the  day  based  on  the  atmospheric 
conditions  and  astronautical  calculations  (Tofsted  1993).  This  model  followed  a 
similar  problem-solving  process;  however,  the  method  used  in  this  iteration  to 
calculate  declination  angle  and  hour  angle  was  slightly  different  in  that  it  used 
equations  from  Jean  Meeus’s  “Astronautical  Algorithms.”  As  such,  an  additional 
0.8  was  subtracted  from  X3,  to  compensate  for  the  difference  in  approach. 

3.  Results 


Results  from  the  model  developed  were  automatically  plotted,  as  part  of  the  SRF 
Model  graphical  user  interface  (GUI).  In  Fig.  2,  the  Surface  Solar  Radiation  Flux 
output  for  a  clear  day  is  displayed. 


ARE2  169  2017-6-18:  Solar  Flux(X  3)W/m  2 


Fig.  2  Model  prediction  2017  June  18 
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These  data  were  visually  compared  to  the  sampled  pyranometer  data  collected  on 
the  clear  day  of  June  18,  2017  (See  Fig.  3).  Figure  3  displays  the  2  pyranometer 
data  time  series  extracted  from  the  Atmospheric  Renewable  Energy  Field  Study  #2 
(ARE2)  measurements.  The  West  Pyranometer  (blue)  was  orientated  to  view  the 
sky  from  a  zenith-facing  angle.  The  East  Pyranometer  (red)  was  aligned  with  a  PV 
Panel’s  south-facing  (32  degrees)  orientation.  The  higher  magnitudes  for  the  angled 
pyranometer  reinforces  the  value  of  angling  the  PV  Panel  for  increased  solar 
radiation  magnitudes. 


ARE2  167  20170618  WSMR,  NM:  Time(hr)  vs.  Solar  Flux(W/m  2) 


Local  Time(hours) 


Fig.  3  Pyranometer  data  170618 

Both  the  SRF  Model  output  and  measured  pyranometer  data  plots  display  a 
Gaussian-shaped  curve  distribution  that  begins  at  roughly  0600  hours  Local  Time 
(LT)  and  ends  at  a  little  past  2000  hours  LT.  The  amplitudes  of  each  rise  up  to 
approximately  970  W/m2  at  just  past  1300  hours  LT.  Note  that  LT  for  the  June  data 
was  Mountain  Daylight  Time  (MDT),  which  explains  the  peak  values  being 
associated  with  1300  LT,  instead  of  the  expected  1200  LT. 

It  can  be  reasonably  asserted  that  the  model  output  aligns  well  with  the  measured 
values  under  clear  atmospheric  conditions.  Moreover,  the  model  appears  to 
coincide  with  the  East  Pyranometer,  which  is  oriented  at  the  same  angle  as  the  solar 
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panel.  This  result  indicates  that  there  is  no  need  to  make  an  adjustment  to  the  model 
to  account  for  the  difference  in  angle. 

4.  Discussion,  Conclusion,  and  Recommendations 

The  initial  SRF  model  based  on  Shapiro’s  model  and  coded/developed  by  David 
Tofsted,  calculated  the  net  solar  radiaiton  flux  and  solar  radiation  flux  without 
scattering,  for  a  single  point  in  time.  For  this  project,  it  was  decided  early  on  that 
the  “updated”  SRF  Model  version  would  be  more  useful  if  the  solar  radiation  flux 
output  was  propagated  throughout  an  entire  solar  day,  and  a  GUI  interface  was 
developed  to  provide  both  a  single  point  output  and  an  all-day  solar  radiation  flux 
plot.  This  strategy  would  more  easily  allow  for  a  comparison  between  model  results 
and  sampled  pyranometer  data,  thus,  preparing  for  the  model  for  its  verification  and 
validation  phase.  When  advancing  from  a  single  time  output  to  an  entire  day,  it 
became  clear  that  the  solar  radiation  flux  without  scattering  did  not  correspond  to 
the  expected  Gaussian  curve.  Since  the  output  did  not  seem  to  relate  to  the  measured 
data  and  the  goal  was  to  employ  relevant  simplicity,  this  feature  was  removed  from 
the  updated  SRF  Model. 

The  model  was  not  adjusted  for  daylight  savings  time,  as  evidenced  by  the  roughly 
1-h  difference  in  solar  noon,  and  a  2-h  difference  in  sunset.  Consequently,  a 
correction  of  60  min  was  included  in  the  calculation  of  solar  time.  Through  testing 
different  dates,  it  was  revealed  that  the  original  declination  angle  and  hour  angle 
calculations  did  not  allow  for  adjustments  of  the  time  of  year  in  the  solar  day.  This 
prompted  declination  angle  and  hour  angle  equation  changes  from  those  algorithms 
used  in  Tofsted’s  code,  to  those  from  Jean  Meeus’s  “Astronautical  Algorithms” 
(see  Section  2).  With  the  equation  changes,  the  model  output  aligned  well  with  the 
available  pyranometer  data.  Unforunately,  most  measured  solar  radiation  data 
available  at  the  conclusion  of  this  task,  were  for  clear  skies.  To  more  strenuously 
test  the  developed  code,  a  variety  of  sky  conditions,  such  as  those  acquired  in  the 
latter  part  of  ARE2,  will  need  to  be  employed. 

For  future  SRF  Model  development,  it  is  recommended  to  develop  a  stepwise 
function  to  facilitate  the  insertion  of  insitu  meteorological  observations  throughout 
a  single  day.  These  perpetual  and  timely  solar  irradiance  updates  could  yield  a  more 
accurate,  dynamic  model.  With  such  a  model,  the  isolated  power  grid  controllers 
would  be  well  poised  for  exploiting  and  optimizing  the  atmospheric  contributions 
to  a  hybrid  power  distribution. 
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Appendix  A.  SolRadCalc  Function 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  [TSpan,  xB,  SolRadO]  =  Sol RadCalc (year,  month,  . 
day,  lat,  long,  tz,  RG,  iRain,  IPU1,  IPUB,  IRT3,  CCV) 

% - 


% - 

%  Function  Title: 

%  Last  Modification: 
%  Code  Author(s) : 

%  Function  Purpose: 

% 

% 


SolRadCalc.m 

170627 

Wal ker/Tofsted 

Compute  the  Net  Downward  Flux  of  Solar  Irradiance 
at  the  Top  of  the  Foliage  Layer  (i.e.  take  account 
of  the  Equivalent  Surface  Reflection) 


%  inputs: 

%  CosZenAng 
%  RG 

%  IRain 

%  IPUl 

% 

%  IPU3 

%  IRT3 

%  CCV 

% 


Current  Cosine  of  the  Solar  Zenith  Angle 
Effective  Reflectivity  of  the  Surface  (Alpha_E) 
Precipitation  Flag  — >  0=No  Precip,  l=Precip. 

Upper  Layer  Clouds  — >  l=Thin,  2=Thick.  (cirrus) 

Middle  Layer  Clouds — >  1, Medium,  always.  (Not  Used) 

Lowest  Layer  Clouds — >  l=Sc/St,  2=Cu/cb.  (IPU3) 

Fog/Smoke  (Ground)  — >  0=None,  l=Present.  (IRT3) 

Fractional  Cloud  Amounts  in  3  Layers  (1,2,3) 
l=Highest  layer,  2=Mid  layer,  3=Low  layer. 


%  Genesis  of  Current  File: 

%  Several  Equations  were  drawn  from... 

%  David  Tofsted,  1993,  "A  Surface  Energy  Budget  Model  Modifying  Heat 

%  Flow  by  Foliage  Effects,"  ARL-TR-60 , 

%  US  Army  Research  Laboratory,  White  Sands  Missile  Range,  NM, 

%  88002-5501. 

%  These  will  be  denoted  by  [TR-60] 


%  Several  Equations  were  drawn  from... 

%  Ralph  Shapiro,  1982,  "Solar  Radiative  Flux  Cal cul cations  From 

%  Standard  Surface  Meteorological  Observations,"  AFGL-TR-82-0039 , 

%  Air  Force  Geophysics  Labratory,  Air  Force  Systems  Command, 

%  United  States  Air  Force,  Hanscom  AFB ,  MA,  01731 

%  These  will  be  denoted  by  [TR-82] 


%  Several  Equations  were  drawn  from... 

%  lean  Meeus,  1999,  "Astronauti cal  Algorithms,"  ISBN  0-943396-61-1 

%  These  will  be  denoted  by  [ISBN  -  1M] 

% - 

%  Calculate  the  cosine  of  the  zenith  angle,  CosZenAng,  the  extraterrestrial 
%  radiaiton  per  unit  horizontal  surface,  X0,  and  the  solar  energy 
%  constant  for  a  quiet  Sun,  SO. 

[X0,  SO,  CosZenAng,  TSpan,  ID]  =  Solarlnfo(year,  month,  day,  lat,  ... 
long,  tz) ; 


% - 

%  Calculate  a  modified  X_3,  that  will  affect  the  solar  radiation  received  at 
%  ground  level  as  a  Cloud  Transmission  Factor  and  the  product  of  the 
%  tranmissivities  of  each  cloud  layer,  Tran_Direct 

[X_3 ,Tran_Di rect]  =  Cl oudTranmulti (CosZenAng,  RG,  IRain,  IPUl,  IPU3,  ... 
IRT3,  CCV); 
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% 

%  Calculate  the  net  solar  flux  at  ground  level,  x3,  and 
%  display  the  solar  radiation  received  at  the  start  of  the  atmosphere, 
%  SolRadO  [TR-60] 
x3  =  x_3 . *x0 ; 

SolRadO  =  S0*ones(l,  250); 


% - 

%  Plot  the  data 

plot(TSpan,  x3,  TSpan,  SolRadO) 

ti  =  sprintf('ARE2  %d  _%d-%d-%d:  Solar  Flux(X_3)  W/mA2 ' ,  JD,  year,  ... 

month,  day); 
title(ti) 

xlabel (' Local  Time(hours) 1 ,  'Fontsize',  18) 
xticks (0:2/24:1) 

xticklabels({'0' , ' 2 ' , ' 4 ' , ' 6 ' , ' 8 ' , ' 10 ' , ' 12 ' , ' 14 ' , ' 16 ' , ' 18 ' , ' 20 ' , ' 22 ' , ' 24 ' }) 
ylabel  ('Solar  Fl ux(w/mA2) ' ,  'Fontsize',  18) 
yticks (0:100: 1500) 
legend('Solar  Flux(Net)',  ... 

'Solar  Constant(Qui et  Sun)',  ... 

'location',  ' northoutside ' ) 
grid  on 
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Appendix  B.  Solarlnfo  Function 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  [XO,  SO,  CosZenAng,  time,  3]  =  Solarlnfo(year,  month,  day,  lat,  ... 
long,  tz) 

% - 


%■ 


%  Function  Title: 

%  Last  Modification: 
%  Code  Author(s) : 

%  Function  Purpose: 

% 

% 

% 


Solarlnfo.m 

170627 

Wal ker/Tofsted 

Compute  the  extraterrestrial  radiation  per  unit 
horizontal  surface,  XO,  and  the  components  used  to 
calculated  it,  SO  and  CosZenAng,  as  shown  in 
Equation  22  of  the  Shapiro  Model. 


%  inputs: 


lat 

local  latitude  in 

degrees 

long 

local  longitude  in 

degrees 

year 

Year (yyyy) 

month 

Month(MM) 

day 

Date(DD) 

tz 

Time  Zone(+  or  -) 

% - 

%Calculate  the  local  Julian  Date  based  on  user  input 
[J,  leap]  =  JulianDate(year,  month,  day); 


% - 

%Calculate  the  number  of  days  since  January  1st,  1900  to  get  to  the  actual 
%  Julian  Day 
datemod  =  year  -  1900; 
date  =  0; 

for  i  =0  :  datemod  -  1 
if  mod(i  ,4)  ==  0 

date  =  date  +  366; 

el  se 

date  =  date  +  365; 

end 

end 

datenum  =  date  +  J ; 

% - 

%Set  up  a  scaled  time  array  for  24  hour  period,  incrementing  by  0.1/24 
time  =  zeros (1,  250); 
time(l)  =  0.1/24; 
for  i  =  2  :  250 

time(i)  =  round(time(i  -  1)  +  (0.1/24),  3); 
time(i  -  1)  =  round(time(i  -1),  3); 

end 


% - 

%Using  Equation  23  of  the  Sharpi ro  model  (d'/d)A2,  a  function  of  the 
%  ell ipti city  of  the  Earth's  orbit  and  the  position  of  the  Earth  in  its 
%  orbit  around  the  Sun,  is  calculated  for  the  current  Julian  Date.  This 
%  variable  is  saved  as  RelDist. 

RelDist  =  1.000140  +  0.016726*cos((2*pi*(J  -  2))/day); 


%■ 
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%Calculate  the  absolute  Julian  Day,  JD,  as  a  function  of  local  Julian  Date 
%  and  the  scaled  time  array.  Convert  Julian  Day  into  Julian  Century,  JC. 

%  The  Geometric  Mean  Longitude  of  the  Sun,  GMLSun ,  and  Geometric  Mean 
%  Anomoly,  GMASun,  are  used  to  calculate  the  position  of  a  body  in  an 
%  elliptical  orbit  as  a  function  of  Julian  Century.  Eccentric  Earth 

%  Orbit,  EEO ,  is  a  measure  of  the  amount  that  Earth  deviates  from  a 

%  perfect  circle  in  its  orbit.  Equation  of  Center  of  the  Sun,  EOC,  is  the 

%  angular  difference  between  the  actual  position  of  a  body  in  its 

%  elliptical  orbit  and  the  position  it  would  occupy  if  its  motion  were 
%  uniform.  Sun's  True  Longitude,  TLSun ,  is  the  ecliptic  longitude  at 
%  which  an  orbiting  body  could  actually  be  found  if  its  inclination  were 

%  zero.  Sun's  Apparent  Longitude,  AppLSun,  is  corrected  for  aberration 

%  and  nutation  as  opposed  to  mean  longitude.  MOblEq  and  oblcorr  account 

%  for  the  angle  at  which  the  Earth  is  tilted,  toward  or  away  from  the  Sun 

%  depending  on  date,  about  its  poles.  [ISBN  -  JM] 
for  i  =  1  :  250 

JD(i)  =  datenum  +  2415018.5  +  time(i)  -  tz/24; 

JC(i)  =  (JD(i)  -  2451545)/36525 ; 

GMLSun (i )  =  mod(280. 46646  +  JC(i)*(36000. 76983  +  JC(i) *0.0003032) , 360) ; 

GMASun (i )  =  357.52911  +  JC(i)*(35999. 05029  -  0.0001537*JC(i)) ; 

EEO(i)  =  0.016708634  -  JC(i)*(0. 000042037  +  0.0000001267*JC(i)) ; 

EOC(i)  =  sin(deg2rad(GMASun(i)))*(l. 914602  -  JC(i)*(0. 004817  +  ... 

0.000014*JC(i)))  +  sin(deg2rad(2*GMASun(i)))*(0. 019993  -  ... 

0.000101*JC(i))  +  sin(deg2rad(3*GMASun(i)))*0. 000289; 

TLSun(i)  =  GMLSun  (i )  +  EOC(i); 

AppLSun(i)  =  TLSun(i)  -  0.00569  -  0.00478*sin(deg2rad(125.04  -  ... 

1934 . 136*JC(i ))) ; 

MOblEq (i)  =  23  +  (26  +  ((21.448  -  JC(i)*(46.815  +  JC(i)*(0. 00059  -  ... 

JC(i )*0. 001813))))/60)/60 ; 

oblcorr(i)  =  MOblEq(i)  +  0.00256*cos(deg2rad(125.04  -  ... 

1934. 136* (JC(i )))) ; 

DecSun(i)  =  rad2deg(asin(sin(deg2rad(oblCorr(i)))*.  .  . 

sin(deg2rad(AppLSun(i))))) ; 

Vary(i)  =  tan(deg2rad(oblCorr(i)/2))A2; 

end 


% - 

%The  Equation  of  Time  is  a  modification  to  the  clock  based  on  the  orbit 
%  of  the  Earth  about  the  Sun  during  a  single  year.  It  modifies  the  time 
%  of  Solar  Noon  based  on  the  Julian  Date.  Sunrise  Angle,  SunRiseAng,  is 

%  the  angle  at  which  the  sun  rises.  Solar  Noon,  SolNoon,  is  when  the  Sun 

%  transits  the  local  celestial  meridian  and  is  at  the  highest  point  in 

%  the  sky.  SunRise  and  sunset  are  the  adjust  times  at  which  the  Sun 

%  crosses  the  Eastern  and  Western  horizons  respectively.  Sunlight 

%  Duration,  SunL,  is  the  total  time  the  Sun  is  above  the  horizon.  True 

%  Solar  time,  SolTime,  is  the  time  based  on  the  apparent  solar  day. 

%  [ISBN  -  JM] 
for  i  =  1  :  250 

EOT(i)  =  4*rad2deg(Vary(i)*sin(2*deg2rad(GMLSun(i)))  -  ... 
2*EEO(i)*sin(deg2rad(GMASun(i)))  +  ... 
4*EEO(i)*Vary(i)*sin(deg2rad(GMASun(i)))*.  .  . 
cos(2*deg2rad(GMLSun(i)))  -  ... 

0. 5*Vary(i)*Vary(i)*sin(4*deg2rad(GMLSun(i)))  -  ... 
1.25*EEO(i)*EEO(i)*si n(2*deg2rad(GMASun(i))))  ; 
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SunRi seAngCi )  =  rad2deg(acos(cos(deg2radC90.833))/(cos(deg2rad(lat))*. . . 
cos(deg2rad(DecSun(i))))  -  tan(deg2rad(lat))*.  .  . 
tan(deg2rad(DecSun(i))))) ; 

SolNoon(i)  =  (720  -  4*long  -  EOT (i )  +  tz*60)/1440; 

SunRise(i)  =  SolNoon(i)  -  4*SunRiseAng(i)/1440; 

SunSet(i)  =  SolNoon(i)  +  4*SunRi  seAng(i)/1440; 

SunL(i )  =  SunRise(i)*8; 

SolTime(i)  =  mod (ti me (i) *1440  +  EOT (i )  +  4*long  -  60*tz,1440); 

end 


%Hour  angle,  HrAng,  is  the  angle  associated  with  the  Sun  at  a  specific  time 
%  and  location.  Negative  angles  here  are  before  solar  noon,  and  positive 
%  angles  after  solar  noon.  Cosine  of  Solar  zenith,  CosZenAng,  angle  based 
%  on  Equation  21  of  the  Shapiro  Model  in  radians.  Additionally,  filter 

%  data  for  times  when  the  zenith  angle  remains  0,  i.e.  when  the  sun  is 

%  not  visible.  [ISBN  -  DM] 
for  i  =  1  :  250 

if  SolTime(i)/4  <  0 

HrAng(i)  =  SolTime(i)/4  +  180; 

el  se 

HrAng(i)  =  SolTime(i)/4  -  180; 

end 

CosZenAng(i)  =  sin(deg2rad(lat))*sin(deg2rad(DecSun(i)))  +  ... 

cos(deg2rad(lat))*cos(deg2rad(DecSun(i)))*cos(deg2rad(HrAng(i))) ; 

end 

% - 

%The  extraterrestrial  radiation  per  unit  horizontal  surface,  X0,  is 
%  calculated  with  Equation  22  of  the  Shapiro  Model,  using  the  solar 

%  constant  SolMax  measured  in  W/mA2,  the  cosine  of  the  solar  zenith  angle 

%  CosZenAng,  and  the  function  of  the  ell ipti city  of  the  Earth's  orbit 
%  and  the  position  of  the  Earth  in  its  orbit  around  the  Sun,  RelDist. 

%  [TR-60] 

SolMax  =  1369.2; 

SO  =  SolMax*RelDistA2 ; 
for  i  =  1  :  250 

if  (CosZenAng (i)  <  0) 

x0(i )  =  0.0; 

el  se 

X0(i)  =  SO. *CosZenAng(i) ; 

end 

end 
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Appendix  C.  CloudTranmulti  Function 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  [X_3,Tran_Di rect]  =  CloudTranmulti (CosZenAng,  RG,  iRain,  IPU1,  ... 
IPU3,  IRT3,  ccv) 


% - 

% - 

%  Function  Title: 

%  Last  Modification: 
%  Code  Author(s) : 

%  Function  Purpose: 

% 


% 


CloudTranmulti  .m 
170627 

Wal ker/Tofsted 

Compute  the  solar  radiation  flux  tranmi  ssivity 
coefficient  based  on  atmospheric  conditions,  cloud 
type,  cloud  density,  and  fractional  cloud  coverage 


%MODEL  TAKEN  FROM: 

%  RALPH  SHAPIRO,  1982,  "SOLAR  RADIATIVE  FLUX  CALCULATIONS  FROM  STANDARD 
%  SURFACE  METEOROLOGICAL  OBSERVATIONS",  AFGL-TR-82-0039 .  AIR  FORCE 
%  GEOPHYSICS  LABORATORY,  HANSCOM  AFB ,  MA  01731 


% 


%  inputs: 

%  CosZenAng 
%  RG 

%  IRain 

%  IPUl 

% 

%  IPU3 

%  IRT  3 

%  CCV 

% 

% - 


Current  Cosine  of  the  Solar  Zenith  Angle 
Effective  Reflectivity  of  the  Surface  (Alpha_E) 
Precipitation  Flag  — >  0=No  Precip,  l=Precip. 

Upper  Layer  Clouds  — >  l=Thin,  2=Thick.  (cirrus) 

Middle  Layer  Clouds — >  1, Medium,  always.  (Not  Used) 

Lowest  Layer  Clouds — >  l=Sc/St,  2=Cu/cb.  (IPU3) 

Fog/Smoke  (Ground)  — >  0=None,  l=Present.  (IRT3) 

Fractional  Cloud  Amounts  in  3  Layers  (1,2,3) 
l=Highest  layer,  2=Mid  layer,  3=Low  layer. 


%Transcribe  the  fractional  cloud 
%  it's  rai  ni  ng(lRai  n  =  1)  then 
%  F(n)  is  the  fractional  cloud 
if  IRain  ==  1 


coverage  from  the  vector  CCV.  However,  if 
all  cloud  layers  are  treated  as  overcast, 
coverage  per  layer. 


Fl  =  1.0; 

F2  =  1.0; 

F3  =  1.0; 

el  se 

Fl  =  CCV(l) 
F2  =  CCV (2) 
F3  =  ccv(3) 


% - 

%Use  diffuse  values  for  layer  2  (midlayer)  when  layer  1  (top  layer) 

%  is  overcast  and  it  is  thick  cirrus,  (obtained  by  setting  ID2  =  1) 
if  (Fl  >=  0.875)  &&  (IPUl  ==  2) 

ID2  =  1; 

el  se 

ID2  =  0; 

end 


% - 

%Use  diffuse  values  for  layer  3  (lowest  layer)  when  either  layer  1 
%  (top  layer)  or  layer  2  (mid  layer)  is  overcast. 

%  (Obtained  by  setting  ID3  =  1). 
if  (ID2  ==  1)  ||  (F2  >  0.875) 
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idB  =  1; 

el  se 

IDB  =  0; 

end 


% - 

%  Reflectivity,  Transmissivity,  and  Weight  Coefficients  (Table  10  -  12) 

% - 

%Tabulate  the  coefficient  matrices,  last  two  values  for  each  variable  in 
%  the  a  and  b  series  are  flipped  to  follow  the  trend  of  the  data  since 

%  there  may  have  been  a  typo  in  the  original  model,  c  series  variable  do 

%  not  change.  [TR-82] 

aO  =  [0.12395,  0.15325,  0.15946,  0.27436,  0.25674,  0.60540,  ... 

0.66152,  0.67072,  0.71214]; 

al  =  [-0.34765,  -0.39620,  -0.42185,  -0.43132,  -0.18077,  -0.5142,  ... 
-0.14863,  -0.13805,  -0.15033]; 

a2  =  [0.39478,  0.42095,  0.48800,  0.26920,  -0.21961,  -0.23389,  ... 

-0.08193,  -0.10895,  0.00696]; 

a3  =  [-0.14627,  -0.14200,  -0.18493,  -0.00447,  0.25272,  0.43648,  ... 

0.13442,  0.09460,  0.03904]; 

bO  =  [0.76977,  0.69318,  0.68679,  0.55336,  0.63547,  0.26498,  ... 

0.19085,  0.17960,  0.13610] ; 

bl  =  [0.49407,  0.68227,  0.71012,  0.61511,  0.35229,  0.66829,  ... 

0.32817,  0.34855,  0.29964]; 

b2  =  [-0.44647,  -0.64289,  -0.71463,  -0.29816,  0.08709,  0.24228,  ... 
-0.08613,  -0.14875,  -0.14041]; 

b3  =  [0.11558,  0.17910,  0.22339,  -0.06663,  -0.22902,  -0.49357,  ... 

-0.08197,  0.01962,  0.00952] ; 
cO  =  [0.575,  1.552,  1.429,  0.858,  2.165]; 
cl  =  [-3.432,  -1.957,  -1.207,  -1.075,  -1.277]; 
c2  =  [1.929,  -1.762,  -2.008,  -0.536,  -3.785]; 
c3  =  [0.842,  2.067,  0.853,  0.750,  2.089]; 
c4  =  [2.693,  0.448,  0.324,  0.322,  -0.387]; 
c5  =  [-1.354,  0.932,  1.582,  0.501,  2.342]; 


% - 

%  Compute  Reflectivity  and  Transmissivity  of  Top  Layer(l).  [TR-82] 

% - 

%r  and  t  for  the  top  layer(l)  are  not  dependant  on  any  other  layers  so  they 
%  will  always  be  calculated  by  the  equations  in  Tables  10  and  11 
for  i  =  1  :  250 

rl(i)  =  a0(l)  +  al(l)*CosZenAng(i)  +  a2(l)*(CosZenAng(i)A2)  +  ... 
a3(l)*(CosZenAng(i)A3) ; 

tl(i)  =  b0(l)  +  bl(l)*CosZenAng(i)  +  b2(l)*(CosZenAng(i)A2)  +  ... 
b3(l)*(CosZenAng(i)A3) ; 

%IPU1  =  1  when  there  are  thin  Ci  clouds  present 
if  IPU1  ==  1 

%Calculate  Reflectivity,  Transmittances,  and  Weighting  Factors  for 
%  overcast  conditions,  i.e.  Fl  >  0.875,  (Table  10-12) 
rhol(i)  =  a0(5)  +  al(5)*CosZenAng(i)  +  a2(5)*(CosZenAng(i)A2)  +  ... 
a3(5)*(CosZenAng(i)A3)  ; 

taul(i)  =  b0(5)  +  bl(5)*CosZenAng(i)  +  b2(5)*(CosZenAng(i)A2)  +  ... 
b3(5)*(CosZenAng(i)A3)  ; 
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Wl(i)  =  cO(l)  +  cl(l)*CosZenAng(i)  +  c2(l)*Fl  +  ... 

c3(l)*Fl*CosZenAng(i)  +  c4(l)*(CosZenAng(i)A2)  +  c5(l)*(FlA2) ; 

%else  IPU  =  2  when  there  are  thick  Ci  clouds  present 
el  se 

%Calculate  Reflectivity,  Transmittances,  and  Weighting  Factors  for 

%  overcast  conditions,  i.e.  Fl  >  0.875,  (Table  10-12) 

rhol(i)  =  a0(6)  +  al(6)*CosZenAng(i)  +  a2(6)*(CosZenAng(i)A2)  +  ... 

aB(6)*(CosZenAng(i)AB) ; 

taul(i)  =  b0(6)  +  bl(6)*CosZenAng(i)  +  b2(6)*(CosZenAng(i)A2)  +  ... 
b3(6)*(CosZenAng(i)A3) ; 

Wl(i)  =  c0(2)  +  cl(2)*CosZenAng(i)  +  c2(2)*Fl  +  ... 

c3(2)*Fl*CosZenAng(i)  +  c4(2)*(CosZenAng(i)A2)  +  c5(2)*(FlA2) ; 

end 

end 

% - 

%  Compute  Reflectivity  and  Transmissivity  of  Middle  Layer(2).  [TR-82] 

% - 

%lmplies  top  layer(l)  is  overcast,  need  to  use  Diffuse  Values  for  middle 
%  layer(2)  found  in  Table  3  of  Shapiro  Model.  Since  middle  layer(2)  has 

%  only  one  cloud  type,  the  same  values  are  used  for  r,  t,  rho,  and  tau. 

%  In  this  case  the  weighting  factor  W  will  determine  overall  reflectivity 
%  and  transmissivity  based  on  the  fractional  cloud  coverage  F(k). 
for  i  =  1  :  250 
if  ID2  ==  1 

r2(i)  =  0.040; 

1 2  C i )  =  0.905; 
rho2(i)  =  0.560; 
tau2(i)  =  0.361; 

W2(i)  =  c0(3)  +  cl(3)*CosZenAng(i)  +  c2(3)*F2  +  ... 

c3(3)*F2*CosZenAng(i)  +  c4(3)*CosZenAng(i)A2  +  c5(3)*F2A2; 

%ID2  =  0  implies  that  the  middle  layer(2)  is  not  diffused 
el  se 

r2(i)  =  a0(2)  +  al(2)*CosZenAng(i)  +  a2(2)*(CosZenAng(i)A2)  +  ... 
a3(2)*(CosZenAng(i)A3) ; 

1 2  C i )  =  b0(2)  +  bl(2)*CosZenAng(i)  +  b2(2)*(CosZenAng(i)A2)  +  ... 
b3(2)*(CosZenAng(i)A3) ; 

%Calculate  Reflectivity,  Transmissivity,  and  Weighting  Factors  for 
%  overcast  conditions,  i.e.  F2  >  0.875,  (Table  10-12) 
rho2(i)  =  a0(7)  +  al(7)*CosZenAng(i)  +  a2(7)*(CosZenAng(i)A2)  +  ... 
a3(7)*(CosZenAng(i)A3)  ; 

tau2(i)  =  b0(7)  +  bl(7)*CosZenAng(i)  +  b2(7)*(CosZenAng(i)A2)  +  ... 
b3(7)*(CosZenAng(i)A3)  ; 

W2(i)  =  c0(3)  +  cl(3)*CosZenAng(i)  +  c2(3)*F2  +  ... 

c3(3)*F2*CosZenAng(i)  +  c4(3)*(CosZenAng(i)A2)  +  c5(3)*(F2A2) ; 

end 

end 


% - 

%  Compute  Reflectivity  and  Transmissivity  of  Bottom  Layer(3).  [TR-82] 

o/o - 
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%ID3  =  1  implies  that  either  the  top  layer(l)  or  the  middle  layer(2)  are 
%  overcast  causing  the  bottom  layer(3)  to  be  diffused 
for  i  =  1  :  250 
if  id3  ==  1 

%IRT3  =  1  implies  that  there  is  fog/smoke  present 
if  IRT3  ==  1 

r3(i)  =  0.116; 
t3(i)  =  0.788; 

%IRT3  =  0  implies  that  there  is  notfog/smoke  present 
el  se 

r3(i)  =  0.045; 
t3(i)  =  0.900; 

end 

%IPU  =  1  implies  that  there  are  Sc/St  clouds  present 
if  IPU3  ==  1 

rho3(i)  =  0.609; 
tau3(i)  =  0.311; 

W3(i)  =  c0(4)  +  cl(4)*CosZenAng(i)  +  c2(4)*F3  +  ... 

c3(4)*F3*CosZenAng(i )  +  c4(4)*(CosZenAng(i )A2)  +  ... 
c5(4)*(f3A2) ; 

%IPU  =  2  implies  that  there  are  Cu/Cb  clouds  present 
el  se 

rho3(i)  =  0.520; 
tau3(i)  =  0.400; 

W3(i)  =  c0(5)  +  cl(5)*CosZenAng(i )  +  c2(5)*F3  +  ... 

c3(5)*F3*CosZenAng(i )  +  c4(5)*(CosZenAng(i )A2)  +  ... 
c5(5)*(F3A2); 

end 

%ID3  =  0  implies  that  the  bottom  layer(3)  is  not  diffused 
el  se 


%IRT3  =  1  implies  that  there  is  fog/smoke  present 
if  IRT3  ==  1 

r3(i)  =  a0(4)  +  al(4)*CosZenAng(i)  +  ... 

a2(4)*(CosZenAng(i)A2)  +  a3(4)*(CosZenAng(i)A3) ; 
t3(i)  =  b0(4)  +  bl(4)*CosZenAng(i)  +  ... 

b2(4)*(CosZenAng(i)A2)  +  b3(4)*(CosZenAng(i)A3) ; 

%IRT3  =  0  implies  that  there  is  not  fog/smoke  present 
el  se 

r3(i)  =  a0(3)  +  al(3)*CosZenAng(i )  +  ... 

a2(3)*(CosZenAng(i)A2)  +  a3(3)*(CosZenAng(i)A3) ; 
t3(i)  =  b0(3)  +  bl(3)*CosZenAng(i )  +  ... 

b2(3)*(CosZenAng(i)A2)  +  b3(3)*(CosZenAng(i)A3) ; 

end 

%IPU  =  1  implies  that  there  are  Sc/St  clouds  present 
if  IPU3  ==  1 

rho3(i)  =  a0(8)  +  al(8)*CosZenAng(i )  +  ... 
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a2(8)*(CosZenAng(i)A2)  +  a3(8)*(CosZenAng(i)A3) ; 

tau3(i)  =  bO(8)  +  bl(8)*CosZenAng(i )  +  ... 

b2(8)*(CosZenAng(i)A2)  +  b3(8)*(CosZenAng(i)A3) ; 

W3(i)  =  c0(4)  +  cl(4)*CosZenAng(i)  +  c2(4)*F3  +  ... 

c3(4)*F3*CosZenAng(i )  +  c4(4)*(CosZenAng(i )A2)  +  ... 
c5(4)*(F3A2) ; 

%IPU  =  2  implies  that  there  are  Cu/Cb  clouds  present 
el  se 

rho3(i)  =  a0(9)  +  al(9)*CosZenAng(i)  +  ... 

a2(9)*(CosZenAng(i)A2)  +  a3(9)*(CosZenAng(i)A3) ; 

tau3(i)  =  b0(9)  +  bl(9)*CosZenAng(i )  +  ... 

b2(9)*(CosZenAng(i)A2)  +  b3(9)*(CosZenAng(i)A3) ; 

W3(i)  =  c0(5)  +  cl(5)*CosZenAng(i )  +  c2(5)*F3  +  ... 

c3(5)*F3*CosZenAngCi )  +  c4(5)*(CosZenAng(i )A2)  +  ... 
c5(5)*(F3A2) ; 

end 

end 

end 


% - 

%  Weighting  Functions,  Eguation  15  [TR-82] 

% - 

%The  weighting  function  value  phi  is  rounded  up  to  1  if  the  fractional 
%  cloud  coverage  is  greater  than  95%,  indicating  that  the  layer  is 
%  severely  overcast.  The  clear  sky  portion  of  the  reflectivity  and 
%  transmissivity  is  ignored  in  this  case, 
for  i  =  1  :  250 
if  (Fl  >  0.95) 

phil(i)  =  1.0; 

el  se 

phil(i)  =  wl(i)*Fl; 

end 

if  (F2  >  0.95) 

phi  2  Ci )  =  1.0; 

el  se 

phi  2 Ci )  =  W2  (i  )*F2  ; 

end 

if  (F3  >  0.95) 
phi  3  Ci )  =  1.0; 

el  se 

phi  3  Ci )  =  w3(i)*F3; 

end 

end 


% - 

%  Reflectivities  and  Tranmissivities  of  the  cloud  layers, 

%  Equations  12  &  13  [TR-82] 

% - 

%Tran_Direct  is  modified  by  reducing  it  by  0.09,  this  accounts  for  the 
%  different  method  used  to  find  CosZenAng  and  HrAng 
for  i  =  1  :  250 
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Rl(i )  =  phil(i) 
Tl(i )  =  phil(i) 
R2(i)  =  phi  2  Ci  D 
T2(i)  =  phi  2(i) 
RB(i)  =  phiB(i) 
T3(i)  =  phi  3(i) 


*rhol(i) 

*taul(i) 

*rho2(i) 

*tau2(i) 

*rho3(i) 

*tau3(i) 


(1 

Cl 

Cl 

Cl 

Cl 

Cl 


philCi))*rlCi) 

PhilCi))*tlCi) 

phi2Ci))*r2Ci) 

Phi2Ci))*t2Ci) 

Phi3Ci))*r3Ci) 

phi3Ci))*t3Ci) 


if  CosZenAngCi)  <=  0 
Tran_Di rectCi)  =  0; 

el  se 

Tran_Di  rectCi)  =  TlCi)*T2Ci)*T3Ci)  ; 

end 


end 


%— 

% 

%— 

for 


Total  Solar  Flux  at  ground  level,  Equation  9  [TR-82] 
i  =  1  :  250 

dlCi)  =  1.0  -  RlCi)"R2Ci); 
d 2 C i )  =  1-0  -  R2Ci)"R3Ci); 
d 3 Ci  )  =  1-0  -  R3Ci)"RG; 

DlCi)  =  dlCi)~d2Ci)  -  RlCi)*R3Ci)*T2Ci)A2; 

D2 Ci )  =  d3Ci)*DlCi)  -  RG*CT3Ci)A2)*CdlCi)"R2Ci)  +  RlCi)*T2Ci)A2)  ; 


x_3Ci)  =  CTran_Di  rectCi)/D2Ci))-  0.08; 


end 
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Intentionally  left  blank. 
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Appendix  D.  JulianDate  Function 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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function  [3D,  leap]  = 

Jul ianDate(yyyy ,  mm,  dd) 

%  Function  Title: 

JulianDate.m 

%  Last  Modification: 

170627 

%  Code  Author(s) : 

Wal ker/Tofsted 

%  Function  Purpose: 

% 

Converts  calander  date  into  Julian  Date  and 
confirms  or  denies  that  it  is  a  leap  year. 

%  inputs: 

%  YYYY  Year 

%  MM  Month 

%  DD  Day 

0/ 

%Error  check  for  date 
if  dd  >  31 

error( ' Entered  invalid  date,  check  user  input') 

end 


% - 

%Set  default  state  to  non-leap  year  to  be  modified  on  a  positive 
%  check  later  in  code 
leap  =  0; 


% - 

%Assign  a  date  number  based  on  the  current  month  entered  with  an 
%  embedded  error  code  to  check  if  a  proper  month  was  entered 
switch  mm 
case  1 

month  =  0; 
case  2 

month  =  0  +  31; 
case  3 

month  =  0  +  31  +  28; 
case  4 

month  =  0  +  31  +  28  +  31; 
case  5 

month  =  0  +  31  +  28  +  31  +  30; 
case  6 

month  =0+31+28+31+30+31; 
case  7 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30; 
case  8 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30  +  31; 
case  9 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30  +  31  +  31; 
case  10 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30  +  31  +  31  +  30; 
case  11 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30  +  31  +  31  +  30  +  31; 
case  12 

month  =  0  +  31  +  28  +  31  +  30  +  31  +  30  +  31  +  31  +  30  +  31  +  30; 
otherwi se 

error( ' Entered  invalid  month,  check  user  input1); 
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end 


% - 

%Julian  dates  are  modified  after  January  in  the  case  of  a  leap  year 
%  if  both  cases  are  positive  here,  the  program  treats  the  date  as 
%  a  leap  year. 

if  (mm  >=  2)  &&  (mod(yyyy,4)  ==  0) 
leap  =  1; 

end 


% - 

%Combine  all  calculated  components  of  the  date  and  return 
JD  =  month  +  leap  +  dd; 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


ARL 

ARE2 

DoD 

EEO 

EOC 

EOT 

ft 

GML 

h 

JC 

JD 

JD 

k 

LT 

MDT 

NWP 

OblCorr 

9 

cpk 

P 

PV 

r 

Rg 

ROTC 

5 

S 


US  Army  Research  Laboratory 

Atmospheric  Renewable  Energy  Field  Study  #2 

Department  of  Defense 

Eccentric  Earth  Orbit 

Equation  of  Center  of  the  Sun 

Equation  of  Time 

fractional  cloud  coverage 

Geometric  Mean  Longitude  of  the  Sun 

Hour  angle 

Julian  Century 

local  Julian  Date 

Julian  Day 

representative  of  the  layer  that  is  being  calculated 
Local  Time 

Mountain  Daylight  Time 

numerical  weather  prediction 

Oblique  correction 

weighting  coefficient  expanded  upon 

weighting  coefficient 

Overcast 

photovoltaic 

clear  layer  reflectivity  coefficient  calculated 
relative  albedo 

Reserved  Officer’s  Training  Corp 
declination  angle 
quiet  Sun  solar  constant 
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SL 

SN 

SR 

SRAng 

SRFM 

SS 

t 

x 

0 

TL 

TST 

W 

X0 

X3 

X3’ 

Y 

Z 


Sunlight  Duration 

Solar  Noon 

Sunrise 

Sunrise  Angle 

Solar  Radiation  Flux  Model 

Sunset 

clear  sky  transmissivity  coefficient  calculated 
overcast  layer  transmissivity  coefficient  calculated 
current  latitude, 

True  Longitude 
True  Solar  time 

bi-quadratic  polynomial  function 
extraterrestrial  radiation 
radiation  flux 

coefficient  of  transmissivity  due  to  atmospheric  conditions 
place  holder  for  the  modification  of  the  oblique  correction 
zenith  angle 
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